Phase separation in the Hubbard model 
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Phase separation in the Hubbard model is investigated with the dynamical cluster approximation. 
We find that it is present in the paramagnetic solution for values of filling smaller than one and at 
finite temperature when a positive next-nearest neighbor hopping is considered. The phase separated 
region is characterized by a mixture of a strongly correlated metallic and Mott insulating phases. 
Our results indicate that phase separation is driven by the formation of doped regions with strong 
antiferromagnetic correlations and low kinetic energy. 



Introduction There is strong experimental evidence 
that high Tc materials are susceptible to charge inhomo- 
geneities, such as stripe&i or checkerboard modulation^. 
This discovery has spurred great theoretical interest in 
phase separation (PS) in models related to the cuprates, 
such as the Hubbard model which is believed to cap- 
ture the low-energy physics of cuprate superconductors. 
It was argued by different authors that the charge in- 
stability displayed as PS in such simple models without 
long-range Coulomb interaction evolves into incommen- 
surate charge ordering when the long-range repulsion is 
considered^. In this paper we present results on PS in the 
Hubbard model. We find that the paramagnetic asym- 
metric Hubbard model near half filling phase separates 
into undoped Mott liquid and doped Mott gas phases. 
The resulting Mott liquid-Mott gas phase diagram bears 
a strong resemblance to that of a classical liquid gas mix- 
ture. 

Phase separation in the Hubbard and in the closely re- 
lated t-J model has been intensively investigated. There 
is a general consensus that a t-J model with a large J/t 
separates into two phases, an undoped antiferromagnet 
(AF) and a hole rich region. However the results for 
realistic J/t < 1 are controversial. Emery et alA, Hell- 
berg et al^ and Gimm et al& report PS for all values 
of J/t. Others authors such as Putikka et alX and Shih 
et al.^ find no PS for small J/t. In the Hubbard model 
with only nearest-neighbor hopping, exact diagonaliza- 
tionS. and Monte Carloi^ calculations show no evidence 
of PS. These numerical results are consistent with the 
analytical results of G. Su's^^, who show that there is no 
phase separation in the particle-hole symmetric Hubbard 
model. However, a large- investigation of this model 
in the infinite U limit shows PS when the next-nearest 
neighbor hopping t' is consideredi^. Phase separation in 
the Hubbard model at small doping was also found in 
a dynamical mean field calculation in the antiferromag- 
netic phasei^ and with variational cluster perturbation 
theoryi^ in the antiferromagnetic and superconducting 
phases. 

Phase separation is believed to be closely related with 
the antiferromagnetic order; a homogeneous doped sys- 
tem is unstable preferring to separate into an undoped 
antiferromagnetic region which lowers the exchange en- 
ergy (maximizes the number of antiferromagnetic bonds) 



and a rich doped phase with low kinetic energy. The 
driving force for PS in a t-J model when J/t is large will 
therefore be the desire to form undoped antiferromag- 
netic regions^. However in the Hubbard model we did not 
find PS for the values of parameters which are optimal 
for antiferromagnetic order in the undoped region. For 
instance, with DCA the maximum Neel temperature in 
the undoped system is obtained for U « 3 /AW, W — 8t 
being the electronic bandwidth, and for t' = 0. The later 
conditions can be understood by noticing that a finite 
t' introduces an antiferromagnetic exchange between the 
same sublattice sites, thus frustrating the antiferromag- 
netism. Nevertheless we find PS only for a U > W and 
a finite next-nearest-neighbor hopping t' ^ 0. Moreover, 
we find PS in the paramagnetic solution which shows that 
short range antiferromagnetic correlations are sufficient 
for the PS to take place. Presumably the PS is driven by 
the formation of weakly doped regions with strong anti- 
ferromagnetic correlations and low kinetic energy. The 
main culprit for the low value of the kinetic energy is the 
parameter t' with the right sign. 

Formalism We use the Dynamical Cluster Approxi- 
mation (DCA)iSii& to explore the possibility of PS in the 
2D Hubbard model, with 



H — Hkin + Hpot 



where 



(1) 

(2) 
(3) 



Here c^jj (creates) destroys an electron with spin a on 
site i and nia- is the corresponding number operator. 
U is the on-site Coulomb repulsion. We consider hop- 
ping t between nearest-neighbors (ij) and hopping t' be- 
tween next- nearest- neighbors {{il)). We show results for 
t — 1, t' = 0.3 and U = 8, which are realistic values 
for cupratesiSiiSiiSi. We find PS for values of the filling 
smaller than one, which for positive <' corresponds to the 
electron doped cuprates. 

The DCA is an extension of the Dynamical Mean Field 
Theory (DMFT)2a. The DMFT maps the lattice prob- 
lem to an impurity embedded self-consistently in a host 
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FIG. 1: Nc = S results. Filling n versus chemical potential 
when T > Tc ~ 0.1 t. Inset: Inverse of the charge suscepti- 
bility Xcharge versus temperature for fixed chemical potential 



and therefore neglects spatial correlations. In the DC A 
we assume that correlations are short-ranged and map 
the original lattice model onto a periodic cluster of size 
Nc = Lc X Lc embedded in a self-consistent host. Thus, 
correlations up to a range ^ ^ Lc are treated accurately, 
while the physics on longer length-scales is described at 
the mean-field level. We solve the cluster problem using 
quantum Monte Carlo (QMC)^^ . The cluster self-energy 
is used to calculate the properties of the host, and this 
procedure is repeated until a self-consistent convergent 
solution is reached. 

Unlike most of the other numerical calculations on PS, 
which study systems with a fixed number of particles, 
our calculations are done in the grand canonical ensem- 
ble and in the thermodynamic limit. Therefore, unlike in 
finite cluster calculations, we do not encounter any par- 
ticular difficulty associated with the small doping regime. 
Phase separation is explored by calculating the filling de- 
pendence on the chemical potential and the charge sus- 
ceptibility (or compressibility), Xcharge = ^■ 

Results First we consider the case of an eight site 
{Nc = 8) cluster. The filling as a function of the chemical 
potential is plotted in Fig.^for different temperatures^. 
Note that that at small doping with lowering temperature 
the charge susceptibility is increasing and diverging at a 
critical point (6c, fic,Tc). The divergence of the charge 
susceptibility is illustrated in the inset. It is a clear in- 
dication that the filling is unstable and the system is 
subject to phase separation into regions with different 
hole density. The critical point is characterized by the 
temperature Tc « 0.10 t and the doping Sc ~ 4.5%. 

For temperatures smaller than Tc and for values of the 
chemical potential close to /ic the DCA calculation pro- 
vides two distinct solutions for the same value of /i. As 
mentioned before, the DCA equations are solved self- 
consistently starting with an initial guess for the self- 
energy, usually zero or that from a larger temperature or 
a perturbation theory result. In most of the situations 



FIG. 2: Nc — S results. Filling n versus chemical potential 
below Tc, at r = 0.071 t. Two solutions describing a hystere- 
sis are found, one incompressible with n ~ 1 (squares) and a 
doped one (circles). Inset: stability of the two solutions ver- 
sus DCA iterations when fj, — 2.96t (middle of the hysteresis, 
corresponding to the dotted line in the main figure). 



an unique solution is obtained independent of the start- 
ing guess. This is the case at doping values far from Sc 
such as 0% (undoped) or 10% doping. However, close 
to fj,c we find that the final solution is dependent on the 
starting point. If one uses as the initial input the self- 
energy corresponding to the undoped solution [n = 1), 
then n versus ^ will look as the upper curve (squares) in 
Fig. 12 On the other hand if the starting self-energy is the 
one corresponding to the large doped solution (n < 1), 
n versus /i will be described by the lower curve (circles) 
in Fig. [3 In both cases, the fully converged self energy 
of the previous point is used to initialize the calculation. 
Thus, below Tc the filling as a function of the chemical 
potential displays a hysteresis. 

Simple thermodynamic ideas may be used to inter- 
pret these results. A hysteresis implies the existence of 
a metastable state and it is observed in many systems 
which suffer a first order transition, a common exam- 
ple being magnetization versus the applied magnetic field 
{M{H)) in magnetic materials. However in the real sys- 
tems, after a sufficient time, the fluctuations always drive 
the system to the stable solution (the equilibrium solu- 
tion) and the hysteresis becomes a discontinuity charac- 
teristic to first order transitions. In our case, due to 
the mean-field coupling of the cluster to the effective 
medium, the hysteresis is stable. This is shown in the 
inset of Fig. |21 where a large number of iterations in the 
self-consistent process is considered. 

By analogy with the liquid-gas system discussed below, 
we label the two states found for T < Tc as Mott liquid 
(ML) and Mott gas (MG). The Mott Hquid is incom- 
pressible and insulating. Both the compressibility and 
doping of the ML are small and decrease with decreasing 
temperature. Its density of states at the Fermi surface 
develops a gap with lowering temperature characteristic 
of an insulator, as seen in Fig.|31-a. The MG is compress- 
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FIG. 3; a) DOS for ML solution at T = 0.077t (dotted line) 
and r = omit (full line), b) DOS for MG solution at T = 
0.077t (dotted line) and T = 0.057^ (full line). 



ible and metallic. The DOS is peaked at the chemical 
potential and increases with lowering temperature, see 
Fig. El-b). Consistent with the narrow peak width, the 
MG is a strongly correlated state with a small value of 
the double occupancy {{rii-^riii) /n « 0.04 at T = 0.077f) 
and strong AF correlations. 

The stable solution below Tc, ML or MG, is the one 
with lower free energy, F = E — jiN — TS. Unfortu- 
nately, due to the mean-field character of the DCA, the 
self-consistent solution is not necessarily the equilibrium 
state, and the QMC method does not allow the calcula- 
tion of the entropy. Therefore the determination of the 
critical ^ where the jump in N{ij) should take place is 
difficult to identify. However, the calculation of the en- 
ergy provides valuable information about the transition 
mechanism. The energy plotted versus jjL displays a cusp 
at /ic when T = Tc (not shown). Below Tc, the energy is 
hysteretic. As can be seen in Fig.01-a at fixed /i the energy 
of the gas phase is much smaller, due to the large gain 
in kinetic energy (see Fig. ^h) produced by the next- 
nearest-neighbor hopping t' as we will discuss. On the 
other hand, the term — /iiV will favor the ML state since 
it has a larger filling. In fact we find that the difference 
between E — fiN for the two solutions is small, with the 
ML state being favored for larger values of /i. When the 
chemical potential is decreased the system will be driven 
to the MG state by both the lower kinetic energy and the 
larger, presumably, entropy characteristic to MG state. 
Therefore, for T < Tc, we expect the jump in n will move 
to lower values of fi as the temperature is lowered. 

One can notice that a phase diagram with these charac- 
teristics bears a striking similarity to the phase diagram 
of a classical liquid gas mixtureSS, where fi plays the role 
of pressure. A cartoon which summarizes our results and 
illustrates this similarity is shown in Fig [S] At high T, 
n versus fj, is linear, since correlations are irrelevant. As 
the temperature is lowered, n{^) becomes nonlinear due 
to correlation effects. At Tc, dn/dfj, diverges. Below Tc 
the hysteresis appears. Upon lowering the temperature 
the hysteresis broadens and the MG (ML) solution shifts 
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FIG. 4: Energy per site versus /i for the two solutions: a) 
total energy E = {H), (see Eq.[TJ at T = 0.077t, b) Kinetic 
energy = {Hki„) (see Eq. EJ at T = 0.077f, c) E - fj.N at 
T = O.Om and T = 0.057t. 
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FIG. 5: Schematic representation of the phase diagram. 



to slightly larger (smaller) dopings. As T — > 0, the en- 
tropy term becomes smaller and the chemical potential 
fic where the jump takes place in the real solution should 
move to smaller values. If a fixed N is imposed when 
T < Tc in the two-phase parameter regime, the system 
will separate into distinct ML and MG regions. 

Our calculations assume a paramagnetic host which im- 
plies that the range of possible AF order is restricted 
to the cluster size. For the Nc = 8 cluster we find PS 
below the AF critical temperature T/v, the temperature 
where the AF spin susceptibility is diverging and the AF 
correlations range reaches the cluster size. Therefore it 
is important to address the role of AF correlations on 
phase separation. For this we investigate the behavior 
of the critical temperatures Tc and T/y when the cluster 
size increases. In the inset of Fig. (Sjone can see that at 
5% doping Tn decreases rapidly with increasing cluster 
size. On the other hand, the Nc = 12 and Nc = 16 site 
clusters display a divergent charge susceptibility roughly 
at the same Tc as the Nc — 8 cluster, Tc « O.lt, as shown 
in Fig. El The rapid decrease of T/v with Nc and the 
fact that Tc is nearly independent on Nc indicates that 
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FIG. 6: Inverse charge susceptibility 1/Xcharge versus temper- 
ature for Nc = S (squares), Nc — 12 (diamonds) and Nc = 16 
(circles) site clusters at 5% doping. Inset: The AF tempera- 
ture Tjv versus cluster size Nc. 



PS may persist in larger clusters at a temperature higher 
than T/v where the range of AF correlations is smaller 
than the cluster size. However, we must mention that 
the calculations on Nc = 12 and Nc = 16 clusters, close 
to the PS temperature, are extremely difficult. This re- 
gion of parameter space is characterized by very strong 
critical behavior (presumably because a larger cluster im- 
plies a weaker hybridization with the effective medium, 
i.e. the results are less "mean-field"), a severe minus sign 
problem, and extremely large auto correlation times be- 
tween measurements. Consequently, the error bar in the 
filling and the charge susceptibility is increasingly large 
for the low temperature points in Fig. and it is dif- 
ficult to obtain converged solutions. Therefore, besides 
a rough estimation of Tc « O.lt it is difficult to make 
other quantitative estimates for the critical parameters. 
Low temperature calculations on larger clusters inside 
the critical region where a hysteresis is expected are not 
possible due to the severe sign problem which appears in 
the QMC calculation. 

We find PS only when the next-nearest-neighbor hop- 
ping t' > and the filling n < 1. A finite t' in Hubbard 
and t-J models is known to give rise to a strong asym- 
metry between electron-doped {f > 0) and hole-doped 
{t' < 0) systemsiSi22i22i2i. In exact diagonalization stud- 
ies on small cluster3i2i2242Si2S42i it was shown that, due 
to the kinetic energy gain, the motion of holes caused 
by a positive t' stabilizes antiferromagnetic configura- 
tions^*'^^. Even though exact diagonalization on systems 
with four holes suggests that t' is not favorable to hole 
clustering-?'—' the tendency to PS was noted in Ref^. 



Our calculations also do not indicate hole clustering but 
rather formation of a 8% — 10% doped state with strong 
AF correlations and low-kinetic energy. Presumably for 
this value of the doping the effect of t' on the kinetic 
energy is the most significant. 

For smaller values of t', PS takes place at lower tem- 
peratures. For instance when t' = O.lt, the system shows 
PS at Tc = 0.055t for the Nc ^ 8 cluster. For U < W we 
found no sign of PS for temperatures above 0.04 t. The 
charge susceptibility behavior suggests that PS is not fa- 
vored when t' < and n < 1, in agreement with exact 
diagonalization results22i2S which show that in this case 
the effect of t' is to push the holes apart from each other. 

Our results imply phase separation into two regions 
with different electronic density. However, even with- 
out considering long range order, we cannot exclude the 
possibility that PS competes with the formation of dif- 
ferent charge patterns such as stripes or checkerboard. 
The investigation of these instabilities would require cal- 
culations on much larger clusters, able to commensurate 
these patterns, which are unfeasible at the moment. 

It would be interesting to investigate the competition 
between PS and d-wave superconductivity in the Hub- 
bard model. However this implies a region of the param- 
eter space not accessible to our method. Calculations on 
clusters larger than Nc = 8 show PS but the sign problem 
precludes access to temperatures where the superconduc- 
tivity is expected. On the other hand, calculations on the 
small 2x2 cluster, where the sign problem is mild, show 
d-wave superconductivity for finite t' but no definite ev- 
idence for PS, even though the charge susceptibility is 
strongly increased when a positive t' is considered^i. 

Conclusions With the DCA we show that the Hub- 
bard model with a positive next-nearest-neighbor hop- 
ping displays PS for values of the filling slightly smaller 
than one. Our results suggest that the PS is driven by the 
desire to form slightly doped (« 8% — 10%) regions with 
low-kinetic energy and strong antiferromagnetic correla- 
tions. The phase diagram is similar to that of the liquid- 
gas mixture, showing a second order critical point and a 
first order transition from a Mott gas to a Mott liquid 
state below Tc- 
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